suppressPackageStartupMessages({
library(Seurat)
library(dplyr)
library(data.table)
library(ggplot2)
library(patchwork)
library(harmony)
library(symphony)
library(presto)
library(singlecellmethods)
library(purrr)
library(pheatmap)
library(ggthemes)
})
set.seed(1)
source("utils.r")
source("myfun.r")
data.path <- "/data/brennerlab/Shani/projects/Treg/data/external_datasets/JIA_Croft/"
list.files(data.path)
metadata <- read.table(paste0(data.path, "object/meta.tsv"))
head(metadata)
# meta.donor <- read.csv(paste0(data.path, "all.meta.CD.csv")) %>% group_by(biosample_id) %>% distinct()
# head(meta.donor)
dim(metadata)
| orig.ident | nCount_RNA | nFeature_RNA | percent.mt | nCount_Protein | nFeature_Protein | T_clonotype_id | T_cdr3s_aa | B_clonotype_id | B_cdr3s_aa | ⋯ | global1 | onset | maxjoint0623 | label15 | simple15 | clusters2312 | oldclusters | label15s | subclusters | main_clusters | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <int> | <int> | <dbl> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | ⋯ | <chr> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | |
| AAACCTGAGACAGGCT-1_1 | PBMC1 | 2437 | 1009 | 2.379975 | 263 | 93 | clonotype3221 | TRB:CAISQTSATYEQYF | NA | NA | ⋯ | T cell | 518 | 1 | CD4+ T | T cells | CD4+ naive/central memory T | CD4+ naive | CD4+ T cells | CD4+ naive/central memory T | CD4+ T cells |
| AAACCTGAGACATAAC-1_1 | PBMC1 | 3180 | 1428 | 7.106918 | 351 | 103 | clonotype7162 | TRB:CASSVDLAGATDTQYF | NA | NA | ⋯ | T cell | 518 | 1 | CD4+ T | T cells | CD4+ KLRB1+ memory T | CD4+ KLRB1+ memory | CD4+ T cells | CD4+ KLRB1+ memory T | CD4+ T cells |
| AAACCTGAGAGGTAGA-1_1 | PBMC1 | 3189 | 1486 | 2.759486 | 289 | 90 | NA | NA | NA | NA | ⋯ | NK/yd T | 518 | 1 | NK cells | NK cells | IFNG+ NK | IFNG+ NK | NK cells | IFNG+ NK | NK cells |
| AAACCTGAGCCCTAAT-1_1 | PBMC1 | 2057 | 884 | 4.472533 | 297 | 89 | clonotype5326 | TRB:CASSYLAGGTNEQFF | NA | NA | ⋯ | T cell | 518 | 1 | CD4+ T | T cells | CD4+ KLRB1+ memory T | CD4+ KLRB1+ memory | CD4+ T cells | CD4+ KLRB1+ memory T | CD4+ T cells |
| AAACCTGAGCCGCCTA-1_1 | PBMC1 | 2645 | 981 | 2.684310 | 292 | 95 | clonotype6751 | TRB:CASSPDRVLSAEAFF | NA | NA | ⋯ | T cell | 518 | 1 | CD4+ T | T cells | CD4+ naive/central memory T | CD4+ naive | CD4+ T cells | CD4+ naive/central memory T | CD4+ T cells |
| AAACCTGAGCGTTCCG-1_1 | PBMC1 | 2306 | 951 | 2.385082 | 223 | 86 | clonotype9493 | TRB:CASSYPRERHSGANVLTF;TRA:CAVNRDDKIIF | NA | NA | ⋯ | T cell | 518 | 1 | Naive CD8+ T | T cells | CD8+ naive T | CD8+ naive | CD8+ T cells | CD8+ naive T | CD8+ T cells |
unique(metadata$main_clusters)
unique(metadata$subclusters)
unique(metadata$label15)
metadata %>% filter(subclusters == "CD4+ FOXP3+ Tregs") %>% dim
metadata %>% filter(main_clusters == "CD4+ T cells") %>% dim
counts <- Read10X(paste0(data.path, "/object/"))
obj <- CreateSeuratObject(counts = counts, project = "JIA_Croft", min.cells=0, min.features=0)
obj
obj@meta.data%>% head
| orig.ident | nCount_RNA | nFeature_RNA | |
|---|---|---|---|
| <fct> | <dbl> | <int> | |
| AAACCTGAGACAGGCT-1_1 | JIA_Croft | 2437 | 1009 |
| AAACCTGAGACATAAC-1_1 | JIA_Croft | 3180 | 1428 |
| AAACCTGAGAGGTAGA-1_1 | JIA_Croft | 3189 | 1486 |
| AAACCTGAGCCCTAAT-1_1 | JIA_Croft | 2057 | 884 |
| AAACCTGAGCCGCCTA-1_1 | JIA_Croft | 2645 | 981 |
| AAACCTGAGCGTTCCG-1_1 | JIA_Croft | 2306 | 951 |
mt <- obj@meta.data %>% tibble::rownames_to_column("cellID") %>%
left_join(metadata %>% tibble::rownames_to_column("cellID"), by= join_by(cellID == cellID))
obj <- AddMetaData(obj, mt)
obj@meta.data %>% head
| orig.ident | nCount_RNA | nFeature_RNA | cellID | orig.ident.x | nCount_RNA.x | nFeature_RNA.x | orig.ident.y | nCount_RNA.y | nFeature_RNA.y | ⋯ | global1 | onset | maxjoint0623 | label15 | simple15 | clusters2312 | oldclusters | label15s | subclusters | main_clusters | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <fct> | <dbl> | <int> | <chr> | <fct> | <dbl> | <int> | <chr> | <int> | <int> | ⋯ | <chr> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | |
| AAACCTGAGACAGGCT-1_1 | JIA_Croft | 2437 | 1009 | AAACCTGAGACAGGCT-1_1 | JIA_Croft | 2437 | 1009 | PBMC1 | 2437 | 1009 | ⋯ | T cell | 518 | 1 | CD4+ T | T cells | CD4+ naive/central memory T | CD4+ naive | CD4+ T cells | CD4+ naive/central memory T | CD4+ T cells |
| AAACCTGAGACATAAC-1_1 | JIA_Croft | 3180 | 1428 | AAACCTGAGACATAAC-1_1 | JIA_Croft | 3180 | 1428 | PBMC1 | 3180 | 1428 | ⋯ | T cell | 518 | 1 | CD4+ T | T cells | CD4+ KLRB1+ memory T | CD4+ KLRB1+ memory | CD4+ T cells | CD4+ KLRB1+ memory T | CD4+ T cells |
| AAACCTGAGAGGTAGA-1_1 | JIA_Croft | 3189 | 1486 | AAACCTGAGAGGTAGA-1_1 | JIA_Croft | 3189 | 1486 | PBMC1 | 3189 | 1486 | ⋯ | NK/yd T | 518 | 1 | NK cells | NK cells | IFNG+ NK | IFNG+ NK | NK cells | IFNG+ NK | NK cells |
| AAACCTGAGCCCTAAT-1_1 | JIA_Croft | 2057 | 884 | AAACCTGAGCCCTAAT-1_1 | JIA_Croft | 2057 | 884 | PBMC1 | 2057 | 884 | ⋯ | T cell | 518 | 1 | CD4+ T | T cells | CD4+ KLRB1+ memory T | CD4+ KLRB1+ memory | CD4+ T cells | CD4+ KLRB1+ memory T | CD4+ T cells |
| AAACCTGAGCCGCCTA-1_1 | JIA_Croft | 2645 | 981 | AAACCTGAGCCGCCTA-1_1 | JIA_Croft | 2645 | 981 | PBMC1 | 2645 | 981 | ⋯ | T cell | 518 | 1 | CD4+ T | T cells | CD4+ naive/central memory T | CD4+ naive | CD4+ T cells | CD4+ naive/central memory T | CD4+ T cells |
| AAACCTGAGCGTTCCG-1_1 | JIA_Croft | 2306 | 951 | AAACCTGAGCGTTCCG-1_1 | JIA_Croft | 2306 | 951 | PBMC1 | 2306 | 951 | ⋯ | T cell | 518 | 1 | Naive CD8+ T | T cells | CD8+ naive T | CD8+ naive | CD8+ T cells | CD8+ naive T | CD8+ T cells |
unique(colnames(obj@meta.data))
unique(obj$global1)
unique(obj$Phase)
unique(obj$TYPE)
unique(obj$sample)
unique(obj$global1)
donor_info <- metadata %>% select(patient, TYPE, sample, sex, batch, age, condition, njoint, ana) %>% unique()
donor_info
| patient | TYPE | sample | sex | batch | age | condition | njoint | ana | |
|---|---|---|---|---|---|---|---|---|---|
| <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <chr> | <int> | <chr> | |
| AAACCTGAGACAGGCT-1_1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N |
| AAACCTGAGACTCGGA-1_2 | 84 | Blood | B84 | F | batch1 | 15.8 | Poly-neg | 24 | N |
| AAACCTGAGACCACGA-1_3 | 95 | Blood | B95 | M | batch1 | 3.7 | P-Oligo | 4 | N |
| AAACCTGAGGTGATTA-1_4 | 811 | Blood | B811 | M | batch2 | 10.8 | P-Oligo | 1 | N |
| AAACCTGAGGATGCGT-1_5 | 88 | Blood | B88 | F | batch2 | 15.7 | RF+Poly | 23 | N |
| AAACCTGAGACAGGCT-1_6 | 98 | Blood | B98 | F | batch2 | 9.3 | P-Oligo | 3 | Y |
| AAACCTGAGGACCACA-1_7 | 97 | Blood | B97 | F | batch2 | 1.2 | Poly-neg | 7 | Y |
| AAACCTGAGAGTGAGA-1_8 | 814 | Blood | B814 | M | batch3 | 14.8 | PsA | 3 | NR |
| AAACCTGAGAGAGCTC-1_9 | 817 | Blood | B817 | F | batch3 | 9.3 | P-Oligo | 1 | Y |
| AAACCTGAGCGTGAAC-1_10 | 91 | SF | SF91 | F | batch1 | 6.9 | P-Oligo | 1 | N |
| AAACCTGAGACACGAC-1_11 | 84 | SF | SF84 | F | batch1 | 15.8 | Poly-neg | 24 | N |
| AAACCTGAGCAACGGT-1_12 | 95 | SF | SF95 | M | batch1 | 3.7 | P-Oligo | 4 | N |
| AAACCTGAGGAATGGA-1_13 | 811 | SF | SF811 | M | batch2 | 10.8 | P-Oligo | 1 | N |
| AAACCTGAGAGCTATA-1_14 | 88 | SF | SF88 | F | batch2 | 15.7 | RF+Poly | 23 | N |
| AAACCTGAGCCACTAT-1_15 | 98 | SF | SF98 | F | batch2 | 9.3 | P-Oligo | 3 | Y |
| AAACCTGAGAATCTCC-1_16 | 97 | SF | SF97 | F | batch2 | 1.2 | Poly-neg | 7 | Y |
| AAACCTGTCAGAGACG-1_17 | 814 | SF | SF814 | M | batch3 | 14.8 | PsA | 3 | NR |
| AAACCTGAGAAGGTTT-1_18 | 817 | SF | SF817 | F | batch3 | 9.3 | P-Oligo | 1 | Y |
| AAACCTGAGAAGGACA-1_19 | 95 | Tissue | T95 | M | batch2 | 3.7 | P-Oligo | 4 | N |
| AAACCTGAGAATTCCC-1_20 | 84 | Tissue | T84 | F | batch2 | 15.8 | Poly-neg | 24 | N |
| AAACCTGAGTTTCCTT-1_21 | 811 | Tissue | T811 | M | batch2 | 10.8 | P-Oligo | 1 | N |
| AAACCTGAGACCTAGG-1_22 | 88 | Tissue | T88 | F | batch2 | 15.7 | RF+Poly | 23 | N |
| AAACCTGAGAAGGTGA-1_23 | 98 | Tissue | T98 | F | batch2 | 9.3 | P-Oligo | 3 | Y |
| AAACCTGAGCAATATG-1_24 | 814 | Tissue | T814 | M | batch3 | 14.8 | PsA | 3 | NR |
| AAACCTGAGAGAGCTC-1_25 | 817 | Tissue | T817 | F | batch3 | 9.3 | P-Oligo | 1 | Y |
| AAACCTGAGACTCGGA-1_26 | 910 | Tissue | T910 | F | batch3 | 3.0 | Ex-Oligo | 1 | Y |
| AAACCTGAGAATTCCC-1_27 | 96 | Tissue | T96 | F | batch3 | 4.8 | PsA | 1 | N |
| AAACCTGAGAGGACGG-1_28 | 99 | Tissue | T99 | F | batch3 | 2.6 | P-Oligo | 2 | NR |
table(donor_info$batch, donor_info$patient)
table(donor_info$patient, donor_info$TYPE)
84 88 91 95 96 97 98 99 811 814 817 910
batch1 2 0 2 2 0 0 0 0 0 0 0 0
batch2 1 3 0 1 0 2 3 0 3 0 0 0
batch3 0 0 0 0 1 0 0 1 0 3 3 1
Blood SF Tissue
84 1 1 1
88 1 1 1
91 1 1 0
95 1 1 1
96 0 0 1
97 1 1 0
98 1 1 1
99 0 0 1
811 1 1 1
814 1 1 1
817 1 1 1
910 0 0 1
# Keep T cells only
Idents(obj) <- "global1"
obj <- subset(obj, ident = "T cell")
rownames(obj)[grep(pattern = "^MT-", rownames(obj))]
obj[["percent.mito"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
VlnPlot(obj, features = "percent.mito", pt.size = 0.0, group.by = "patient") + theme(legend.position = "none") + ggtitle(NULL) +
xlab(NULL) + ylab("% mitochondrial genes") + geom_hline(yintercept=15)
VlnPlot(obj, features = "percent.mito", pt.size = 0.0, group.by = "batch") + theme(legend.position = "none") + ggtitle(NULL) +
xlab(NULL) + ylab("% mitochondrial genes") + geom_hline(yintercept=15)
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
#Gene number
VlnPlot(obj, features = "nFeature_RNA", pt.size = 0.0, group.by = "patient") + theme(legend.position = "none") + ggtitle(NULL) +
xlab(NULL) + ylab("Number of Genes")
VlnPlot(obj, features = "nFeature_RNA", pt.size = 0.0, group.by = "batch") + theme(legend.position = "none") + ggtitle(NULL) +
xlab(NULL) + ylab("Number of Genes")
VlnPlot(obj, features = "nCount_RNA", pt.size = 0.0, group.by = "patient") + theme(legend.position = "none") + ggtitle(NULL) +
xlab(NULL) + ylab("nUMI")
VlnPlot(obj, features = "nCount_RNA", pt.size = 0.0, group.by = "batch") + theme(legend.position = "none") + ggtitle(NULL) +
xlab(NULL) + ylab("nUMI")
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
fig.size(8,8)
suppressMessages(
obj@meta.data %>%
ggplot(aes(x=nCount_RNA, y=nFeature_RNA, color=percent.mito)) +
geom_point() +
scale_colour_gradient(low = "gray90", high = "black") +
facet_wrap(~orig.ident) +
stat_smooth(method=lm) +
scale_x_log10() +
scale_y_log10() +
theme_classic() +
geom_vline(xintercept = 500) +
geom_hline(yintercept = 250) +
theme(text = element_text(size = 14))
)
#Data already filtered
# obj <- subset(obj, subset = nCount_RNA > 500 & nFeature_RNA > 250 & percent.mito < 20)
#initial dimensionality reduction on all samples to identify contaminating clusters
obj <- NormalizeData(obj, normalization.method = "LogNormalize", scale.factor = 10000, verbose = F) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>%
ScaleData() %>%
RunPCA(verbose = F)
Centering and scaling data matrix
set.seed(1)
obj <- RunHarmony(obj, group.by.vars=c("batch", "patient", "TYPE"), reduction.use = "pca",
plot_convergence = TRUE, max_iter = 10,
early_stop = F)
Transposing data matrix Initializing state using k-means centroids initialization Harmony 1/10 Harmony 2/10 Harmony 3/10 Harmony 4/10 Harmony 5/10 Harmony 6/10 Harmony 7/10 Harmony 8/10 Harmony 9/10 Harmony 10/10
obj
An object of class Seurat 30803 features across 100098 samples within 1 assay Active assay: RNA (30803 features, 2000 variable features) 3 layers present: counts, data, scale.data 2 dimensional reductions calculated: pca, harmony
set.seed(1)
obj <- Run_uwot_umap(obj)
obj <- FindClusters(obj, graph.name = 'humap_fgraph', resolution = 1, verbose = TRUE)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 100098 Number of edges: 1234828 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.8233 Number of communities: 14 Elapsed time: 14 seconds
fig.size(8,8)
DimPlot(obj, label = T, repel = T, pt.size = 0.5, label.size = 6) + #scale_color_tableau("Tableau 20") +
theme(text = element_text(size = 20))
DimPlot(obj, group.by = "batch")
DimPlot(obj, group.by = "patient")
DimPlot(obj, group.by = "TYPE")
DimPlot(obj, group.by = "condition")
fig.size(8,16)
DimPlot(obj, group.by = "subclusters")
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE` Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
Rasterizing points since number of points exceeds 100,000. To disable this behavior set `raster=FALSE`
fig.size(12, 20)
FeaturePlot(obj, c("CD3E", "CD4", "FOXP3", "IL2RA", "CTLA4", "IKZF2", "CD8A", "CXCR6"), ncol = 4, raster=T)
FeaturePlot(obj, c("CD3E", "CD4", "FOXP3", "IL2RA", "CTLA4", "IKZF2", "CD8A", "CXCR6"), ncol = 4, order = T, raster=T)
fig.size(4,10)
VlnPlot(obj, "CD4", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "CD8A", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "CD3E", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "CD3D", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "FOXP3", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "IL2RA", group.by = "seurat_clusters", raster=T)
VlnPlot(obj, "CXCL13", group.by = "seurat_clusters", raster=T)
table(obj$subclusters, obj$seurat_clusters)
0 1 2 3 4 5 6 7 8 9
Activated NK-like T 107 7 9 10 1034 10 7632 1 1 6
CD4+ FOXP3+ Tregs 3 38 97 12 3 37 4 113 4327 4
CD4+ KLRB1+ memory T 303 8761 654 71 40 158 42 527 69 76
CD4+ naive/central memory T 30 3670 9539 578 253 9681 34 112 32 8
CD56+bright NK 1 4 1 0 0 1 8 0 0 0
CD8+ GZMB+/GZMK+ memory T 7870 6 1 26 1 1 62 12 3 1548
CD8+ GZMK+ memory T 5853 125 186 105 756 12 286 21 3 215
CD8+ naive T 113 40 708 109 8401 278 354 4 0 13
CXCL13+ TpH 72 1774 333 161 11 157 18 5327 1209 449
Cycling Myeloid 0 0 1 0 0 0 8 0 0 0
Cycling T 4 7 0 1 2 0 50 1 1 3
IFNG+ NK 72 0 0 0 0 0 5 0 0 3
MAIT cells 3 33 4 3 0 0 0 0 0 0
NK-like ILC1 12 0 14 79 1 0 1 0 0 0
THEMIS+ IL7R+ ILC 235 73 635 9945 67 197 62 46 82 30
Vdelta2 yd 27 4 0 0 1 0 10 0 0 0
10 11 12 13
Activated NK-like T 14 0 0 19
CD4+ FOXP3+ Tregs 1 80 13 6
CD4+ KLRB1+ memory T 0 0 131 40
CD4+ naive/central memory T 3 7 136 64
CD56+bright NK 2 0 0 9
CD8+ GZMB+/GZMK+ memory T 4 0 10 16
CD8+ GZMK+ memory T 24 0 53 23
CD8+ naive T 7 0 57 22
CXCL13+ TpH 2 29 90 17
Cycling Myeloid 0 0 0 0
Cycling T 0 0 0 0
IFNG+ NK 5 0 0 34
MAIT cells 0 0 0 0
NK-like ILC1 1183 0 0 66
THEMIS+ IL7R+ ILC 367 990 4 112
Vdelta2 yd 0 0 0 5
#subset
tregs <- subset(obj, ident = c(8))
tregs
An object of class Seurat 30803 features across 5727 samples within 1 assay Active assay: RNA (30803 features, 2000 variable features) 3 layers present: counts, data, scale.data 3 dimensional reductions calculated: pca, harmony, humap
tregs <- NormalizeData(tregs, normalization.method = "LogNormalize", scale.factor = 10000, verbose = F) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>%
ScaleData() %>%
RunPCA(verbose = F)
Centering and scaling data matrix
set.seed(1)
tregs <- RunHarmony(tregs, group.by.vars=c("batch", "patient", "TYPE"), reduction.use = "pca",
plot_convergence = TRUE, max_iter = 10,
early_stop = T)
Transposing data matrix Initializing state using k-means centroids initialization Harmony 1/10 Harmony 2/10 Harmony 3/10 Harmony 4/10 Harmony converged after 4 iterations
set.seed(1)
tregs <- Run_uwot_umap(tregs)
tregs <- FindClusters(tregs, graph.name = 'humap_fgraph', resolution = 0.8, verbose = TRUE)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 5727 Number of edges: 69517 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.6827 Number of communities: 9 Elapsed time: 0 seconds
fig.size(8,8)
DimPlot(tregs, label = T, repel = T, pt.size = 0.5, label.size = 6) +
scale_color_tableau("Tableau 20") + theme(text = element_text(size = 20))
DimPlot(tregs, group.by = "batch")
DimPlot(tregs, group.by = "patient")
DimPlot(tregs, group.by = "TYPE")
DimPlot(tregs, group.by = "condition")
fig.size(8,16)
DimPlot(tregs, group.by = "subclusters")
fig.size(12, 15)
FeaturePlot(tregs, c("CD3E", "CD4", "FOXP3", "IL2RA", "CTLA4", "IKZF2", "AREG", "TCF7", "CXCR6"), ncol = 3)
fig.size(4,10)
VlnPlot(tregs, "FOXP3", group.by = "seurat_clusters")
VlnPlot(tregs, "IL2RA", group.by = "seurat_clusters")
VlnPlot(tregs, "AREG", group.by = "seurat_clusters")
VlnPlot(tregs, "TCF7", group.by = "seurat_clusters")
# comapre to source annotations
table(tregs$subclusters, tregs$seurat_clusters)
0 1 2 3 4 5 6 7 8
Activated NK-like T 0 0 0 0 0 0 1 0 0
CD4+ FOXP3+ Tregs 629 903 872 576 802 362 123 35 25
CD4+ KLRB1+ memory T 27 18 5 12 0 1 6 0 0
CD4+ naive/central memory T 28 3 0 1 0 0 0 0 0
CD8+ GZMB+/GZMK+ memory T 0 0 0 0 0 0 3 0 0
CD8+ GZMK+ memory T 0 0 0 1 0 1 1 0 0
CXCL13+ TpH 455 213 99 283 8 103 41 5 2
Cycling T 0 0 0 0 1 0 0 0 0
THEMIS+ IL7R+ ILC 11 0 1 0 2 66 2 0 0
Tregs.markers <- wilcoxauc(tregs, "seurat_clusters")
options(repr.matrix.max.cols=150, repr.matrix.max.rows=200)
Tregs.markers %>%
group_by(group) %>%
filter(padj < 0.05 & auc > 0.6) %>%
arrange(group, desc(logFC)) %>%
#slice_max(n = 40, order_by = logFC)%>%
dplyr::select(feature, group) %>%
mutate(row = row_number()) %>%
tidyr::pivot_wider(names_from = group, values_from = feature) -> Treg.show
Treg.show[1:50,]
| row | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|---|
| <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> |
| 1 | LEF1 | KLF2 | LGALS1 | FOS | TNFRSF18 | MTRNR2L12 | CCL5 | NR4A1 | TRBV7-6 |
| 2 | RPS12 | EMP3 | HLA-DRB1 | JUNB | TNFRSF4 | NEAT1 | NKG7 | DUSP2 | TRBV7-4 |
| 3 | EEF1B2 | TXNIP | S100A4 | DUSP1 | TIGIT | SYNE2 | CD8A | CD69 | TRBV7-2 |
| 4 | TCF7 | TLE5 | CD74 | CREM | LINC01943 | MT-ATP6 | CTSW | FOS | TRBV7-9 |
| 5 | RPS3A | CD52 | S100A11 | ZNF331 | CTLA4 | MALAT1 | CD8B | NFKBID | TRAV38-2DV8 |
| 6 | RPS13 | LIME1 | LGALS3 | FOSB | BATF | MT-CO1 | GZMA | TNFRSF4 | KLRB1 |
| 7 | RPL32 | SELPLG | S100A6 | DUSP2 | CTSC | ARHGAP15 | GZMK | ZFP36L1 | GZMA |
| 8 | RPS8 | S100A10 | VIM | ZFP36 | CARD16 | FOXP1 | CST7 | CTLA4 | STX10 |
| 9 | RPS5 | S100A6 | CRIP1 | DNAJB1 | IL2RA | RNF213 | LAG3 | ZFP36 | PIK3CG |
| 10 | RPS18 | SH3BGRL3 | S100A10 | TNFAIP3 | FOXP3 | MT-ND5 | GNLY | PIM3 | CLSTN3 |
| 11 | RPL13 | TMSB10 | ANXA2 | RGCC | CD74 | MT-CO3 | KLRK1 | RGCC | FBXL14 |
| 12 | RPL30 | B2M | HLA-DRA | CXCR4 | SPOCK2 | SKAP1 | CCL4 | EGR3 | PRKAB1 |
| 13 | TXNIP | HLA-C | CLIC1 | NR4A2 | SRGN | CASK | GZMB | EGR1 | NA |
| 14 | RPL22 | RPL14 | ANXA1 | JUND | IL2RB | SMCHD1 | ZFP36L2 | NR4A3 | NA |
| 15 | RPL3 | RPL18A | UCP2 | CD69 | TYMP | LINC02694 | ANXA1 | PHLDA1 | NA |
| 16 | RPS27 | NA | HLA-DPA1 | JUN | GAPDH | MT-ND6 | CD96 | BATF | NA |
| 17 | RPL39 | NA | MYL6 | DNAJA1 | TNFRSF1B | MBNL1 | THEMIS | TNF | NA |
| 18 | MAL | NA | TIMP1 | LMNA | ENO1 | MT-CYB | CD7 | SH2D2A | NA |
| 19 | RPS6 | NA | EMP3 | YPEL5 | DUSP4 | MT-CO2 | GZMM | ICOS | NA |
| 20 | CCR7 | NA | CD99 | SLC2A3 | CXCR6 | PHACTR2 | PDE3B | KLF6 | NA |
| 21 | RPL21 | NA | GZMA | TENT5C | PKM | CHST11 | AOAH | BTG2 | NA |
| 22 | SELL | NA | HLA-DPB1 | RGS1 | LAYN | EPB41 | PARP8 | NFKBIA | NA |
| 23 | RPS23 | NA | LY6E | PPP1R15A | SAT1 | MT-ND1 | HCST | SRGN | NA |
| 24 | RPS21 | NA | SH3BGRL3 | BTG2 | CST7 | LPP | PTMS | CSRNP1 | NA |
| 25 | RPL34 | NA | HLA-DRB5 | DUSP4 | TPI1 | MT-ND4L | GZMH | GEM | NA |
| 26 | RPLP2 | NA | GAPDH | SELENOK | PHACTR2 | USP15 | STK17A | BIRC3 | NA |
| 27 | RPL18A | NA | CD2 | KLF6 | CD2 | RABGAP1L | S100A6 | TNFRSF18 | NA |
| 28 | GIMAP7 | NA | TXN | PHLDA1 | ACP5 | AKAP13 | SP100 | JUND | NA |
| 29 | RPL5 | NA | CTSC | FTH1 | CD7 | GLCCI1 | LYAR | IER2 | NA |
| 30 | RPS29 | NA | TYMP | TSC22D3 | IGFLR1 | CAMK4 | FXYD2 | CREM | NA |
| 31 | RPL11 | NA | CKLF | HSPH1 | PHLDA1 | PDE3B | ADGRE5 | NR4A2 | NA |
| 32 | RPL9 | NA | ARPC1B | SRGN | MAGEH1 | DPYD | CD99 | TNFRSF1B | NA |
| 33 | RPL10A | NA | COTL1 | UBE2S | PRDM1 | GPRIN3 | CBLB | DNAJA1 | NA |
| 34 | RPS25 | NA | IL2RG | HSPA8 | CACYBP | LDLRAD4 | CD81 | JUNB | NA |
| 35 | RPL35A | NA | ACTG1 | NR4A1 | AC017002.3 | TNIK | PRF1 | PTPN7 | NA |
| 36 | RPS27A | NA | ENO1 | HSPA1B | CD27 | IQGAP2 | TCF25 | GNG2 | NA |
| 37 | RPL29 | NA | CORO1A | HSP90AA1 | IL2RG | MT-ND3 | PIK3R5 | REL | NA |
| 38 | RPS14 | NA | PLP2 | SOCS1 | DNPH1 | PIK3R5 | PITPNC1 | TBC1D4 | NA |
| 39 | RPL36 | NA | MYL12A | BTG1 | PGAM1 | PACS1 | RUNX3 | ARID5B | NA |
| 40 | RPL37 | NA | AHNAK | SAMSN1 | STK17B | KLF12 | LITAF | NAMPT | NA |
| 41 | RPL19 | NA | TPM3 | TAGAP | GBP5 | LRBA | APOBEC3G | RNF19A | NA |
| 42 | KLF2 | NA | HCST | H3F3B | LGALS3 | PRKCH | CYBA | HSP90AA1 | NA |
| 43 | RPL36A | NA | CD52 | HSP90AB1 | ICOS | MT-ND4 | ITGAE | RILPL2 | NA |
| 44 | RPS3 | NA | PKM | CCNH | MYL6 | MT-ND2 | ID2 | PPP1R15A | NA |
| 45 | RPL10 | NA | ACTB | NR4A3 | CBLB | SAMHD1 | CTSD | NEDD9 | NA |
| 46 | RPL38 | NA | CXCR3 | MYADM | IL12RB2 | NA | CXCR6 | NFKB1 | NA |
| 47 | RPL8 | NA | HLA-DQB1 | PTPN22 | RNF213 | NA | TLN1 | ARF6 | NA |
| 48 | EEF1G | NA | CYBA | IRF1 | GBP2 | NA | MBP | MTRNR2L12 | NA |
| 49 | RPL12 | NA | LSP1 | HSPE1 | PTPRJ | NA | EMB | STAT3 | NA |
| 50 | EEF1A1 | NA | RNF213 | NAMPT | LAG3 | NA | GSTP1 | RAB8B | NA |
table(tregs$seurat_clusters)
0 1 2 3 4 5 6 7 8 9 10 11 1538 1224 1178 1102 1018 970 576 506 434 233 152 130
consistent with the strategy for other datasets
# 6 is CD8
# 8 is contamination TRBV genes
tregs <- subset(tregs, ident = c(6,8), invert = T)
#initial dimensionality reduction on all samples to identify contaminating clusters
tregs <- NormalizeData(tregs, normalization.method = "LogNormalize", scale.factor = 10000, verbose = F) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>%
ScaleData() %>%
RunPCA(verbose = F)
Centering and scaling data matrix
tregs
An object of class Seurat 30803 features across 5523 samples within 1 assay Active assay: RNA (30803 features, 2000 variable features) 3 layers present: counts, data, scale.data 3 dimensional reductions calculated: pca, harmony, humap
tregs@meta.data %>% head
| orig.ident | nCount_RNA | nFeature_RNA | cellID | orig.ident.x | nCount_RNA.x | nFeature_RNA.x | orig.ident.y | nCount_RNA.y | nFeature_RNA.y | percent.mt | nCount_Protein | nFeature_Protein | T_clonotype_id | T_cdr3s_aa | B_clonotype_id | B_cdr3s_aa | scrublet_score | integrated_snn_res.0.1 | integrated_snn_res.0.15 | seurat_clusters | integrated_snn_res.0.12 | percent.ribro | percent.Hb | S.Score | G2M.Score | Phase | old.ident | patient | TYPE | sample | sex | batch | age | condition | njoint | ana | global1 | onset | maxjoint0623 | label15 | simple15 | clusters2312 | oldclusters | label15s | subclusters | main_clusters | percent.mito | humap_fgraph_res.1 | humap_fgraph_res.0.8 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <fct> | <dbl> | <int> | <chr> | <fct> | <dbl> | <int> | <chr> | <int> | <int> | <dbl> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <int> | <int> | <fct> | <int> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | <chr> | <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <chr> | <int> | <chr> | <chr> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <dbl> | <fct> | <fct> | |
| AAAGATGAGCAGCGTA-1_1 | JIA_Croft | 1350 | 760 | AAAGATGAGCAGCGTA-1_1 | JIA_Croft | 1350 | 760 | PBMC1 | 1350 | 760 | 3.1111111 | 340 | 99 | clonotype4891 | TRB:CASSLAIRNYGYTF | NA | NA | 0.06086957 | 1 | 1 | 0 | 1 | 0.1507320 | 0 | 0.01850598 | 0.068792887 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ KLRB1+ memory T | CD4+ KLRB1+ memory | CD4+ T cells | CD4+ KLRB1+ memory T | CD4+ T cells | 3.1111111 | 8 | 0 |
| AAAGATGTCCTAGTGA-1_1 | JIA_Croft | 1485 | 773 | AAAGATGTCCTAGTGA-1_1 | JIA_Croft | 1485 | 773 | PBMC1 | 1485 | 773 | 1.4814815 | 292 | 93 | clonotype77 | TRA:CAVRLSWGKLQF | NA | NA | 0.08163265 | 1 | 1 | 1 | 1 | 0.1495413 | 0 | -0.04059452 | 0.008258175 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 1.4814815 | 8 | 1 |
| AAAGTAGCACTAGTAC-1_1 | JIA_Croft | 1567 | 847 | AAAGTAGCACTAGTAC-1_1 | JIA_Croft | 1567 | 847 | PBMC1 | 1567 | 847 | 2.1059349 | 330 | 100 | clonotype1715 | TRB:CASSSQNRAEQYF | NA | NA | 0.05263158 | 1 | 1 | 0 | 1 | 0.1298440 | 0 | 0.01361041 | 0.082733395 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CXCL13+ TpH | CXCL13+ TpH | CD4+ T cells | CXCL13+ TpH | CD4+ T cells | 2.1059349 | 8 | 0 |
| AACCATGGTCGCATCG-1_1 | JIA_Croft | 2589 | 1048 | AACCATGGTCGCATCG-1_1 | JIA_Croft | 2589 | 1048 | PBMC1 | 2589 | 1048 | 0.9269988 | 278 | 89 | clonotype751 | TRB:CSAEETGPETQYF;TRA:CAVTGGTSYGKLTF | NA | NA | 0.15625000 | 1 | 1 | 0 | 1 | 0.1489572 | 0 | 0.04661555 | -0.045203442 | S | S | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 0.9269988 | 8 | 0 |
| AACCATGGTTGCTCCT-1_1 | JIA_Croft | 2049 | 1102 | AACCATGGTTGCTCCT-1_1 | JIA_Croft | 2049 | 1102 | PBMC1 | 2049 | 1102 | 1.2201074 | 279 | 90 | clonotype9315 | TRB:CASRSTGGGWSNPYEQYF;TRA:CVVNPGRGGAQKLVF | NA | NA | 0.03914591 | 1 | 1 | 1 | 1 | 0.1104492 | 0 | -0.04041617 | -0.052987582 | G1 | G1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 1.2201074 | 8 | 1 |
| AACGTTGAGGATGGAA-1_1 | JIA_Croft | 4757 | 1769 | AACGTTGAGGATGGAA-1_1 | JIA_Croft | 4757 | 1769 | PBMC1 | 4757 | 1769 | 1.8499054 | 293 | 92 | clonotype1286 | TRB:CASSPILGGPTEAFF;TRB:CSARDPRQGYEQFF;TRA:CAASHNNNDMRF;TRA:CAVMEYGNKLVF | NA | NA | 0.09497207 | 1 | 1 | 0 | 1 | 0.1185773 | 0 | -0.02773367 | -0.003626403 | G1 | G1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CXCL13+ TpH | CXCL13+ TpH | CD4+ T cells | CXCL13+ TpH | CD4+ T cells | 1.8499054 | 8 | 0 |
fig.size(4,4)
set.seed(1)
tregs <- RunHarmony(tregs, group.by.vars=c("batch", "patient", "TYPE"), reduction.use = "pca",
plot_convergence = TRUE, max_iter = 10,
early_stop = T)
# ElbowPlot(obj, ndims = 50, reduction = "harmony")
# ElbowPlot(obj, ndims = 50, reduction = "pca")
Transposing data matrix Initializing state using k-means centroids initialization Harmony 1/10 Harmony 2/10 Harmony 3/10 Harmony 4/10 Harmony converged after 4 iterations
set.seed(1)
tregs <- Run_uwot_umap(tregs)
tregs <- FindClusters(tregs, graph.name = 'humap_fgraph', resolution = 0.8, verbose = TRUE)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 5523 Number of edges: 67522 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.6758 Number of communities: 8 Elapsed time: 0 seconds
fig.size(6,6)
DimPlot(tregs, label = T, repel = T, pt.size = 0.5, label.size = 6) +
scale_color_tableau("Tableau 20") + theme(text = element_text(size = 20))
DimPlot(tregs, group.by = "batch")
DimPlot(tregs, group.by = "patient")
DimPlot(tregs, group.by = "TYPE")
DimPlot(tregs, group.by = "condition")
fig.size(8,16)
DimPlot(tregs, group.by = "subclusters")
fig.size(12, 12)
FeaturePlot(tregs, c("FOXP3", "IL2RA", "CTLA4", "IKZF2", "CXCR6", "AREG", "MKI67", "TCF7", "CXCL13"), ncol = 3)
FeaturePlot(tregs, order = T, c("FOXP3", "IL2RA", "CTLA4", "IKZF2", "CXCR6", "AREG", "MKI67", "TCF7", "CXCL13"), ncol = 3)
fig.size(5, 30)
FeaturePlot(tregs, order = T, c("FOXP3"), split.by = "seurat_clusters", ncol = 4)
FeaturePlot(tregs, order = F, c("FOXP3"), split.by = "seurat_clusters", ncol = 4)
fig.size(2,10)
VlnPlot(tregs, "FOXP3", group.by = "seurat_clusters")
VlnPlot(tregs, "IL2RA", group.by = "seurat_clusters")
VlnPlot(tregs, "CD4", group.by = "seurat_clusters")
VlnPlot(tregs, "CD8A", group.by = "seurat_clusters")
VlnPlot(tregs, "CXCL13", group.by = "seurat_clusters")
saveRDS(tregs, "/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs.RDS")
# treg <- readRDS("/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/CD_Tregs.RDS")
ref <- "/data/brennerlab/Shani/projects/Treg/analysis/integrated/integrated_Tregs_SymphonyRef.rds"
tregRef<- readRDS(ref)
query = mapQuery(tregs@assays$RNA@counts, # query gene expression (genes x cells)
tregs@meta.data, # query metadata (cells x attributes)
tregRef, # Symphony reference object
vars = c("batch", "patient", "TYPE"), #"PubID", "Site", "Type" # Query batch variables to harmonize over (NULL treats query as one batch)
do_normalize = TRUE, # perform log(CP10k) normalization on query (set to FALSE if already normalized)
do_umap = TRUE) # project query cells into reference UMAP
Normalizing Scaling and synchronizing query gene expression Found 1940 out of 2000 reference variable genes in query dataset Project query cells using reference gene loadings Clustering query cells to reference centroids Correcting query batch effects UMAP All done!
summary(query)
Length Class Mode exp 170124969 dgCMatrix S4 meta_data 54 data.frame list Z 276150 -none- numeric Zq_pca 276150 -none- numeric R 552300 -none- numeric Xq 44184 -none- numeric umap 11046 -none- numeric
fig.size(4,6)
qp <- query$umap %>% as.data.frame %>% cbind(query$meta_data) %>%
ggplot(aes(x=UMAP1, y=UMAP2, color=seurat_clusters)) + geom_point(size = 0.2) +
theme_bw(base_size = 18)
qp
fig.size(4,10)
qp + facet_wrap(facets = vars(TYPE))
qp + facet_wrap(facets = vars(patient))
# Predict query cell types
query <- knnPredict(query, tregRef,
train_labels = tregRef$meta_data$cell.states,
k = 30,
confidence = TRUE) # calculate prediction confidence
query$meta_data %>% head
| orig.ident | nCount_RNA | nFeature_RNA | cellID | orig.ident.x | nCount_RNA.x | nFeature_RNA.x | orig.ident.y | nCount_RNA.y | nFeature_RNA.y | percent.mt | nCount_Protein | nFeature_Protein | T_clonotype_id | T_cdr3s_aa | B_clonotype_id | B_cdr3s_aa | scrublet_score | integrated_snn_res.0.1 | integrated_snn_res.0.15 | seurat_clusters | integrated_snn_res.0.12 | percent.ribro | percent.Hb | S.Score | G2M.Score | Phase | old.ident | patient | TYPE | sample | sex | batch | age | condition | njoint | ana | global1 | onset | maxjoint0623 | label15 | simple15 | clusters2312 | oldclusters | label15s | subclusters | main_clusters | percent.mito | humap_fgraph_res.1 | humap_fgraph_res.0.8 | cell.states | cell.states.knn.prob | donorID | tissue | cell_type_pred_knn | cell_type_pred_knn_prob | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <fct> | <dbl> | <int> | <chr> | <fct> | <dbl> | <int> | <chr> | <int> | <int> | <dbl> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <int> | <int> | <fct> | <int> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | <chr> | <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <chr> | <int> | <chr> | <chr> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <dbl> | <fct> | <fct> | <fct> | <dbl> | <int> | <chr> | <fct> | <dbl> | |
| AAAGATGAGCAGCGTA-1_1 | JIA_Croft | 1350 | 760 | AAAGATGAGCAGCGTA-1_1 | JIA_Croft | 1350 | 760 | PBMC1 | 1350 | 760 | 3.1111111 | 340 | 99 | clonotype4891 | TRB:CASSLAIRNYGYTF | NA | NA | 0.06086957 | 1 | 1 | 0 | 1 | 0.1507320 | 0 | 0.01850598 | 0.068792887 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ KLRB1+ memory T | CD4+ KLRB1+ memory | CD4+ T cells | CD4+ KLRB1+ memory T | CD4+ T cells | 3.1111111 | 8 | 0 | Naive Treg | 0.4333333 | 91 | Blood | Naive Treg | 0.4333333 |
| AAAGATGTCCTAGTGA-1_1 | JIA_Croft | 1485 | 773 | AAAGATGTCCTAGTGA-1_1 | JIA_Croft | 1485 | 773 | PBMC1 | 1485 | 773 | 1.4814815 | 292 | 93 | clonotype77 | TRA:CAVRLSWGKLQF | NA | NA | 0.08163265 | 1 | 1 | 0 | 1 | 0.1495413 | 0 | -0.04059452 | 0.008258175 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 1.4814815 | 8 | 0 | Naive Treg | 0.5666667 | 91 | Blood | CD25int Treg | 0.5666667 |
| AAAGTAGCACTAGTAC-1_1 | JIA_Croft | 1567 | 847 | AAAGTAGCACTAGTAC-1_1 | JIA_Croft | 1567 | 847 | PBMC1 | 1567 | 847 | 2.1059349 | 330 | 100 | clonotype1715 | TRB:CASSSQNRAEQYF | NA | NA | 0.05263158 | 1 | 1 | 1 | 1 | 0.1298440 | 0 | 0.01361041 | 0.082733395 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CXCL13+ TpH | CXCL13+ TpH | CD4+ T cells | CXCL13+ TpH | CD4+ T cells | 2.1059349 | 8 | 1 | Naive Treg | 0.6000000 | 91 | Blood | Naive Treg | 0.6333333 |
| AACCATGGTCGCATCG-1_1 | JIA_Croft | 2589 | 1048 | AACCATGGTCGCATCG-1_1 | JIA_Croft | 2589 | 1048 | PBMC1 | 2589 | 1048 | 0.9269988 | 278 | 89 | clonotype751 | TRB:CSAEETGPETQYF;TRA:CAVTGGTSYGKLTF | NA | NA | 0.15625000 | 1 | 1 | 1 | 1 | 0.1489572 | 0 | 0.04661555 | -0.045203442 | S | S | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 0.9269988 | 8 | 1 | Naive Treg | 0.9333333 | 91 | Blood | Naive Treg | 0.9333333 |
| AACCATGGTTGCTCCT-1_1 | JIA_Croft | 2049 | 1102 | AACCATGGTTGCTCCT-1_1 | JIA_Croft | 2049 | 1102 | PBMC1 | 2049 | 1102 | 1.2201074 | 279 | 90 | clonotype9315 | TRB:CASRSTGGGWSNPYEQYF;TRA:CVVNPGRGGAQKLVF | NA | NA | 0.03914591 | 1 | 1 | 6 | 1 | 0.1104492 | 0 | -0.04041617 | -0.052987582 | G1 | G1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 1.2201074 | 8 | 6 | CD25int Treg | 0.5666667 | 91 | Blood | ISG Treg | 0.5333333 |
| AACGTTGAGGATGGAA-1_1 | JIA_Croft | 4757 | 1769 | AACGTTGAGGATGGAA-1_1 | JIA_Croft | 4757 | 1769 | PBMC1 | 4757 | 1769 | 1.8499054 | 293 | 92 | clonotype1286 | TRB:CASSPILGGPTEAFF;TRB:CSARDPRQGYEQFF;TRA:CAASHNNNDMRF;TRA:CAVMEYGNKLVF | NA | NA | 0.09497207 | 1 | 1 | 1 | 1 | 0.1185773 | 0 | -0.02773367 | -0.003626403 | G1 | G1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CXCL13+ TpH | CXCL13+ TpH | CD4+ T cells | CXCL13+ TpH | CD4+ T cells | 1.8499054 | 8 | 1 | Naive Treg | 0.7666667 | 91 | Blood | Naive Treg | 0.7333333 |
fig.size(8, 12)
qp <- query$umap %>% as.data.frame %>% cbind(query$meta_data) %>%
ggplot(aes(x=UMAP1, y=UMAP2, color=cell_type_pred_knn)) + geom_point(size = 0.8) +
theme_bw(base_size = 18) + guides(colour = guide_legend(override.aes = list(size=2))) +
scale_color_manual(values = cell.state.colors$cell.states) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +coord_fixed()
qp
fig.size(10,18)
qp + facet_wrap(facets = vars(cell_type_pred_knn))
fig.size(6, 16)
qp + facet_wrap(facets = vars(TYPE))
qp + facet_wrap(facets = vars(batch))
qp + facet_wrap(facets = vars(patient))
qp + facet_wrap(facets = vars(age))
query$meta_data %>% select(subclusters, cell_type_pred_knn, cell_type_pred_knn_prob) %>% arrange(cell_type_pred_knn_prob) %>% head
table(query$meta_data$subclusters, query$meta_data$cell_type_pred_knn)
| subclusters | cell_type_pred_knn | cell_type_pred_knn_prob | |
|---|---|---|---|
| <chr> | <fct> | <dbl> | |
| GTGCTTCCATCGGACC-1_12 | CD4+ KLRB1+ memory T | AREG+ Treg | 0.2333333 |
| ACATACGCAAGCTGAG-1_14 | CXCL13+ TpH | CD25int Treg | 0.2333333 |
| TTCGAAGAGCGAAGGG-1_14 | CD4+ FOXP3+ Tregs | CD25highCXCR6+ Treg | 0.2580645 |
| TCAGCTCTCGTCTGAA-1_10 | CXCL13+ TpH | CD161+mem. Treg | 0.2666667 |
| CCAGCGATCTGGCGAC-1_11 | CXCL13+ TpH | CD25high Treg | 0.2666667 |
| CATATGGCAGTTTACG-1_14 | CD4+ FOXP3+ Tregs | AREG+ Treg | 0.2666667 |
Naive Treg CD25int Treg CD25high Treg
CD4+ FOXP3+ Tregs 285 516 1557
CD4+ KLRB1+ memory T 8 18 11
CD4+ naive/central memory T 15 8 0
CD8+ GZMK+ memory T 0 0 0
CXCL13+ TpH 98 271 226
Cycling T 0 0 0
THEMIS+ IL7R+ ILC 16 11 17
CD25highCXCR6+ Treg AREG+ Treg TNFAIP3+ Treg
CD4+ FOXP3+ Tregs 990 212 358
CD4+ KLRB1+ memory T 1 15 2
CD4+ naive/central memory T 0 8 0
CD8+ GZMK+ memory T 0 1 1
CXCL13+ TpH 59 332 94
Cycling T 0 0 0
THEMIS+ IL7R+ ILC 20 11 0
CD161+mem. Treg ISG Treg GZM Treg Prolif.
CD4+ FOXP3+ Tregs 61 179 14 7
CD4+ KLRB1+ memory T 5 0 3 0
CD4+ naive/central memory T 1 0 0 0
CD8+ GZMK+ memory T 0 0 0 0
CXCL13+ TpH 41 38 6 1
Cycling T 0 0 0 1
THEMIS+ IL7R+ ILC 2 2 0 1
fig.size(8, 12)
qp + facet_wrap(facets = vars(subclusters, TYPE))
# qp + facet_wrap(facets = vars(orig.ident, Layer))
mt<- tregs@meta.data %>% tibble::rownames_to_column("rownames") %>%
left_join(query$meta_data %>% tibble::rownames_to_column("rownames"), by = join_by(rownames))
mt %>% head
| rownames | orig.ident.x.x | nCount_RNA.x.x | nFeature_RNA.x.x | cellID.x | orig.ident.x.x.x | nCount_RNA.x.x.x | nFeature_RNA.x.x.x | orig.ident.y.x | nCount_RNA.y.x | nFeature_RNA.y.x | percent.mt.x | nCount_Protein.x | nFeature_Protein.x | T_clonotype_id.x | T_cdr3s_aa.x | B_clonotype_id.x | B_cdr3s_aa.x | scrublet_score.x | integrated_snn_res.0.1.x | integrated_snn_res.0.15.x | seurat_clusters.x | integrated_snn_res.0.12.x | percent.ribro.x | percent.Hb.x | S.Score.x | G2M.Score.x | Phase.x | old.ident.x | patient.x | TYPE.x | sample.x | sex.x | batch.x | age.x | condition.x | njoint.x | ana.x | global1.x | onset.x | maxjoint0623.x | label15.x | simple15.x | clusters2312.x | oldclusters.x | label15s.x | subclusters.x | main_clusters.x | percent.mito.x | humap_fgraph_res.1.x | humap_fgraph_res.0.8.x | cell.states.x | cell.states.knn.prob.x | donorID.x | tissue.x | orig.ident.y.y | nCount_RNA.y.y | nFeature_RNA.y.y | cellID.y | orig.ident.x.y | nCount_RNA.x.y | nFeature_RNA.x.y | orig.ident.y.y.y | nCount_RNA.y.y.y | nFeature_RNA.y.y.y | percent.mt.y | nCount_Protein.y | nFeature_Protein.y | T_clonotype_id.y | T_cdr3s_aa.y | B_clonotype_id.y | B_cdr3s_aa.y | scrublet_score.y | integrated_snn_res.0.1.y | integrated_snn_res.0.15.y | seurat_clusters.y | integrated_snn_res.0.12.y | percent.ribro.y | percent.Hb.y | S.Score.y | G2M.Score.y | Phase.y | old.ident.y | patient.y | TYPE.y | sample.y | sex.y | batch.y | age.y | condition.y | njoint.y | ana.y | global1.y | onset.y | maxjoint0623.y | label15.y | simple15.y | clusters2312.y | oldclusters.y | label15s.y | subclusters.y | main_clusters.y | percent.mito.y | humap_fgraph_res.1.y | humap_fgraph_res.0.8.y | cell.states.y | cell.states.knn.prob.y | donorID.y | tissue.y | cell_type_pred_knn | cell_type_pred_knn_prob | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <fct> | <dbl> | <int> | <chr> | <fct> | <dbl> | <int> | <chr> | <int> | <int> | <dbl> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <int> | <int> | <fct> | <int> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | <chr> | <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <chr> | <int> | <chr> | <chr> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <dbl> | <fct> | <fct> | <fct> | <dbl> | <int> | <chr> | <fct> | <dbl> | <int> | <chr> | <fct> | <dbl> | <int> | <chr> | <int> | <int> | <dbl> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <int> | <int> | <fct> | <int> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | <chr> | <int> | <chr> | <chr> | <chr> | <chr> | <dbl> | <chr> | <int> | <chr> | <chr> | <int> | <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <dbl> | <fct> | <fct> | <fct> | <dbl> | <int> | <chr> | <fct> | <dbl> | |
| 1 | AAAGATGAGCAGCGTA-1_1 | JIA_Croft | 1350 | 760 | AAAGATGAGCAGCGTA-1_1 | JIA_Croft | 1350 | 760 | PBMC1 | 1350 | 760 | 3.1111111 | 340 | 99 | clonotype4891 | TRB:CASSLAIRNYGYTF | NA | NA | 0.06086957 | 1 | 1 | 0 | 1 | 0.1507320 | 0 | 0.01850598 | 0.068792887 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ KLRB1+ memory T | CD4+ KLRB1+ memory | CD4+ T cells | CD4+ KLRB1+ memory T | CD4+ T cells | 3.1111111 | 8 | 0 | Naive Treg | 0.4333333 | 91 | Blood | JIA_Croft | 1350 | 760 | AAAGATGAGCAGCGTA-1_1 | JIA_Croft | 1350 | 760 | PBMC1 | 1350 | 760 | 3.1111111 | 340 | 99 | clonotype4891 | TRB:CASSLAIRNYGYTF | NA | NA | 0.06086957 | 1 | 1 | 0 | 1 | 0.1507320 | 0 | 0.01850598 | 0.068792887 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ KLRB1+ memory T | CD4+ KLRB1+ memory | CD4+ T cells | CD4+ KLRB1+ memory T | CD4+ T cells | 3.1111111 | 8 | 0 | Naive Treg | 0.4333333 | 91 | Blood | Naive Treg | 0.4333333 |
| 2 | AAAGATGTCCTAGTGA-1_1 | JIA_Croft | 1485 | 773 | AAAGATGTCCTAGTGA-1_1 | JIA_Croft | 1485 | 773 | PBMC1 | 1485 | 773 | 1.4814815 | 292 | 93 | clonotype77 | TRA:CAVRLSWGKLQF | NA | NA | 0.08163265 | 1 | 1 | 0 | 1 | 0.1495413 | 0 | -0.04059452 | 0.008258175 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 1.4814815 | 8 | 0 | Naive Treg | 0.5666667 | 91 | Blood | JIA_Croft | 1485 | 773 | AAAGATGTCCTAGTGA-1_1 | JIA_Croft | 1485 | 773 | PBMC1 | 1485 | 773 | 1.4814815 | 292 | 93 | clonotype77 | TRA:CAVRLSWGKLQF | NA | NA | 0.08163265 | 1 | 1 | 0 | 1 | 0.1495413 | 0 | -0.04059452 | 0.008258175 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 1.4814815 | 8 | 0 | Naive Treg | 0.5666667 | 91 | Blood | CD25int Treg | 0.5666667 |
| 3 | AAAGTAGCACTAGTAC-1_1 | JIA_Croft | 1567 | 847 | AAAGTAGCACTAGTAC-1_1 | JIA_Croft | 1567 | 847 | PBMC1 | 1567 | 847 | 2.1059349 | 330 | 100 | clonotype1715 | TRB:CASSSQNRAEQYF | NA | NA | 0.05263158 | 1 | 1 | 1 | 1 | 0.1298440 | 0 | 0.01361041 | 0.082733395 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CXCL13+ TpH | CXCL13+ TpH | CD4+ T cells | CXCL13+ TpH | CD4+ T cells | 2.1059349 | 8 | 1 | Naive Treg | 0.6000000 | 91 | Blood | JIA_Croft | 1567 | 847 | AAAGTAGCACTAGTAC-1_1 | JIA_Croft | 1567 | 847 | PBMC1 | 1567 | 847 | 2.1059349 | 330 | 100 | clonotype1715 | TRB:CASSSQNRAEQYF | NA | NA | 0.05263158 | 1 | 1 | 1 | 1 | 0.1298440 | 0 | 0.01361041 | 0.082733395 | G2M | G2M | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CXCL13+ TpH | CXCL13+ TpH | CD4+ T cells | CXCL13+ TpH | CD4+ T cells | 2.1059349 | 8 | 1 | Naive Treg | 0.6000000 | 91 | Blood | Naive Treg | 0.6333333 |
| 4 | AACCATGGTCGCATCG-1_1 | JIA_Croft | 2589 | 1048 | AACCATGGTCGCATCG-1_1 | JIA_Croft | 2589 | 1048 | PBMC1 | 2589 | 1048 | 0.9269988 | 278 | 89 | clonotype751 | TRB:CSAEETGPETQYF;TRA:CAVTGGTSYGKLTF | NA | NA | 0.15625000 | 1 | 1 | 1 | 1 | 0.1489572 | 0 | 0.04661555 | -0.045203442 | S | S | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 0.9269988 | 8 | 1 | Naive Treg | 0.9333333 | 91 | Blood | JIA_Croft | 2589 | 1048 | AACCATGGTCGCATCG-1_1 | JIA_Croft | 2589 | 1048 | PBMC1 | 2589 | 1048 | 0.9269988 | 278 | 89 | clonotype751 | TRB:CSAEETGPETQYF;TRA:CAVTGGTSYGKLTF | NA | NA | 0.15625000 | 1 | 1 | 1 | 1 | 0.1489572 | 0 | 0.04661555 | -0.045203442 | S | S | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 0.9269988 | 8 | 1 | Naive Treg | 0.9333333 | 91 | Blood | Naive Treg | 0.9333333 |
| 5 | AACCATGGTTGCTCCT-1_1 | JIA_Croft | 2049 | 1102 | AACCATGGTTGCTCCT-1_1 | JIA_Croft | 2049 | 1102 | PBMC1 | 2049 | 1102 | 1.2201074 | 279 | 90 | clonotype9315 | TRB:CASRSTGGGWSNPYEQYF;TRA:CVVNPGRGGAQKLVF | NA | NA | 0.03914591 | 1 | 1 | 6 | 1 | 0.1104492 | 0 | -0.04041617 | -0.052987582 | G1 | G1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 1.2201074 | 8 | 6 | CD25int Treg | 0.5666667 | 91 | Blood | JIA_Croft | 2049 | 1102 | AACCATGGTTGCTCCT-1_1 | JIA_Croft | 2049 | 1102 | PBMC1 | 2049 | 1102 | 1.2201074 | 279 | 90 | clonotype9315 | TRB:CASRSTGGGWSNPYEQYF;TRA:CVVNPGRGGAQKLVF | NA | NA | 0.03914591 | 1 | 1 | 6 | 1 | 0.1104492 | 0 | -0.04041617 | -0.052987582 | G1 | G1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CD4+ FOXP3+ Tregs | CD4+ FOXP3+ Tregs | CD4+ T cells | CD4+ FOXP3+ Tregs | CD4+ T cells | 1.2201074 | 8 | 6 | CD25int Treg | 0.5666667 | 91 | Blood | ISG Treg | 0.5333333 |
| 6 | AACGTTGAGGATGGAA-1_1 | JIA_Croft | 4757 | 1769 | AACGTTGAGGATGGAA-1_1 | JIA_Croft | 4757 | 1769 | PBMC1 | 4757 | 1769 | 1.8499054 | 293 | 92 | clonotype1286 | TRB:CASSPILGGPTEAFF;TRB:CSARDPRQGYEQFF;TRA:CAASHNNNDMRF;TRA:CAVMEYGNKLVF | NA | NA | 0.09497207 | 1 | 1 | 1 | 1 | 0.1185773 | 0 | -0.02773367 | -0.003626403 | G1 | G1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CXCL13+ TpH | CXCL13+ TpH | CD4+ T cells | CXCL13+ TpH | CD4+ T cells | 1.8499054 | 8 | 1 | Naive Treg | 0.7666667 | 91 | Blood | JIA_Croft | 4757 | 1769 | AACGTTGAGGATGGAA-1_1 | JIA_Croft | 4757 | 1769 | PBMC1 | 4757 | 1769 | 1.8499054 | 293 | 92 | clonotype1286 | TRB:CASSPILGGPTEAFF;TRB:CSARDPRQGYEQFF;TRA:CAASHNNNDMRF;TRA:CAVMEYGNKLVF | NA | NA | 0.09497207 | 1 | 1 | 1 | 1 | 0.1185773 | 0 | -0.02773367 | -0.003626403 | G1 | G1 | 91 | Blood | B91 | F | batch1 | 6.9 | P-Oligo | 1 | N | T cell | 518 | 1 | CD4+ T | T cells | CXCL13+ TpH | CXCL13+ TpH | CD4+ T cells | CXCL13+ TpH | CD4+ T cells | 1.8499054 | 8 | 1 | Naive Treg | 0.7666667 | 91 | Blood | Naive Treg | 0.7333333 |
## Add umap and cell state predictions to seurat obj.
tregs@reductions[["umap"]] <- CreateDimReducObject(embeddings = query$umap[rownames(tregs@meta.data),], key = "UMAP_", assay = DefaultAssay(tregs))
tregs <- AddMetaData(tregs, mt$cell_type_pred_knn, "cell.states")
tregs <- AddMetaData(tregs, mt$cell_type_pred_knn_prob, "cell.states.knn.prob")
# saveRDS(query, "/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs_symphony.RDS")
# saveRDS(tregs, "/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs_symphony_Seurat.RDS")
tregs <- readRDS("/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs_symphony_Seurat.RDS")
query <- readRDS( "/data/brennerlab/Shani/projects/Treg/analysis/Treg_nonRA/JIA_Croft_Tregs_symphony.RDS")
tregs
An object of class Seurat 30803 features across 5523 samples within 1 assay Active assay: RNA (30803 features, 2000 variable features) 3 layers present: counts, data, scale.data 4 dimensional reductions calculated: pca, harmony, humap, umap
query$meta_data %>% names()
fig.size(6,6)
count.groups <- query$meta_data %>% group_by(patient) %>% mutate(totalP = n()) %>% ungroup()%>%
select(cell_type_pred_knn, TYPE, patient, totalP) %>%
group_by(cell_type_pred_knn, TYPE, patient) %>% mutate(nType = n()) %>%
mutate(propType = nType/ totalP) %>% ungroup() %>%
group_by(cell_type_pred_knn, TYPE) %>% summarize(prop.average = mean(propType))
count.groups
ggplot(count.groups, aes(x = TYPE, y = prop.average, fill = cell_type_pred_knn)) +
geom_col(position = "fill") + scale_fill_manual(values = cell.state.colors.v2$cell.states) + theme_base(base_size = 20)
`summarise()` has grouped output by 'cell_type_pred_knn'. You can override using the `.groups` argument.
| cell_type_pred_knn | TYPE | prop.average |
|---|---|---|
| <fct> | <chr> | <dbl> |
| Naive Treg | Blood | 0.098614599 |
| Naive Treg | SF | 0.009955783 |
| Naive Treg | Tissue | 0.016194000 |
| CD25int Treg | Blood | 0.084877277 |
| CD25int Treg | SF | 0.098925497 |
| CD25int Treg | Tissue | 0.112756653 |
| CD25high Treg | Blood | 0.049476911 |
| CD25high Treg | SF | 0.290763208 |
| CD25high Treg | Tissue | 0.200442745 |
| CD25highCXCR6+ Treg | Blood | 0.002112659 |
| CD25highCXCR6+ Treg | SF | 0.191614706 |
| CD25highCXCR6+ Treg | Tissue | 0.158955875 |
| AREG+ Treg | Blood | 0.005190914 |
| AREG+ Treg | SF | 0.023143855 |
| AREG+ Treg | Tissue | 0.205168262 |
| TNFAIP3+ Treg | SF | 0.006382039 |
| TNFAIP3+ Treg | Tissue | 0.135917068 |
| CD161+mem. Treg | Blood | 0.008174522 |
| CD161+mem. Treg | SF | 0.015064627 |
| CD161+mem. Treg | Tissue | 0.017886946 |
| ISG Treg | Blood | 0.006153404 |
| ISG Treg | SF | 0.070800637 |
| ISG Treg | Tissue | 0.007396221 |
| GZM Treg | Blood | 0.004656989 |
| GZM Treg | SF | 0.004756350 |
| GZM Treg | Tissue | 0.007300054 |
| Prolif. | Blood | 0.001792115 |
| Prolif. | SF | 0.002680202 |
| Prolif. | Tissue | 0.005463594 |
fig.size(6,6)
count.groups <- query$meta_data %>% select(cell_type_pred_knn, condition) %>% group_by(cell_type_pred_knn, condition) %>% count()
fig.size(6,6)
count.groups <- query$meta_data %>% group_by(patient) %>% mutate(totalP = n()) %>% ungroup()%>%
select(cell_type_pred_knn, condition, patient, totalP) %>%
group_by(cell_type_pred_knn, condition, patient) %>% mutate(nType = n()) %>%
mutate(propType = nType/ totalP) %>% ungroup() %>%
group_by(cell_type_pred_knn, condition) %>% summarize(prop.average = mean(propType))
ggplot(count.groups, aes(x = condition, y = prop.average, fill = cell_type_pred_knn)) +
geom_col(position = "fill") + scale_fill_manual(values = cell.state.colors.v2$cell.states) +
theme_base(base_size = 20) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplot(count.groups, aes(x = condition, y = prop.average, fill = cell_type_pred_knn)) +
geom_col(position = "stack") + scale_fill_manual(values = cell.state.colors.v2$cell.states) +
theme_base(base_size = 20) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
`summarise()` has grouped output by 'cell_type_pred_knn'. You can override using the `.groups` argument.
table(query$meta_data$condition, query$meta_data$patient)
84 88 91 95 96 97 98 99 811 814 817 910
Ex-Oligo 0 0 0 0 0 0 0 0 0 0 0 552
P-Oligo 0 0 558 592 0 0 820 270 281 0 496 0
Poly-neg 563 0 0 0 0 263 0 0 0 0 0 0
PsA 0 0 0 0 194 0 0 0 0 238 0 0
RF+Poly 0 696 0 0 0 0 0 0 0 0 0 0
table(query$meta_data$cell_type_pred_knn, query$meta_data$patient)
84 88 91 95 96 97 98 99 811 814 817 910
Naive Treg 50 38 73 90 4 42 54 2 9 33 16 11
CD25int Treg 139 88 100 146 22 29 68 22 31 48 33 98
CD25high Treg 183 207 200 171 96 59 337 74 185 44 168 87
CD25highCXCR6+ Treg 23 195 102 96 40 85 202 57 36 45 156 33
AREG+ Treg 84 74 5 56 16 7 24 61 2 24 31 195
TNFAIP3+ Treg 41 66 1 7 4 3 68 47 2 34 64 118
CD161+mem. Treg 4 11 12 11 8 8 18 3 9 5 16 5
ISG Treg 37 16 60 9 2 30 45 1 6 5 7 1
GZM Treg 1 0 4 6 0 0 2 2 1 0 3 4
Prolif. 1 1 1 0 2 0 2 1 0 0 2 0
obj$condition %>% unique()
fig.size(6,6)
query$meta_data <- query$meta_data %>% group_by(patient, T_clonotype_id) %>% mutate(cloneSize = n()) %>%
mutate(cloneSize = if_else(is.na(T_clonotype_id), 0, cloneSize))%>%
ungroup()
qp <- query$umap %>% as.data.frame %>% cbind(query$meta_data) %>% arrange(cloneSize) %>%
ggplot(aes(x=UMAP1, y=UMAP2, color=as.factor(cloneSize))) + geom_point(size = 0.1, position = ) +
theme_bw(base_size = 18) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +coord_fixed()
qp
metadata$cloneSize %>% table
. 0 1 2 3 4 5 6 7 8 9 12 23 26 340 3603 822 297 136 100 48 49 40 27 12 23 26
fig.size(6,8)
metadata <- query$meta_data %>%
mutate(cloneCategory = ifelse(is.na(T_clonotype_id), "No TCR information",
ifelse(cloneSize == "1", "Singlet",
ifelse(cloneSize %in% c(2:3), "Small (1 < X <= 3)",
ifelse(cloneSize %in% c(4:10), "Medium (3 < X <= 10)",
ifelse(cloneSize %in% c(11:30), "Large (10 < X <= 30)",
"Hyperexpanded (30 < X)")))))) %>%
mutate(clonal = ifelse(cloneSize > 1, "Clonal", "Not Clonal"))
# metadata$cloneCategory <- tidyr::replace_na(metadata$cloneCategory, "No TCR information")
metadata$cloneCategory <- factor(metadata$cloneCategory, levels = c("Hyperexpanded (30 < X)",
"Large (10 < X <= 30)",
"Medium (3 < X <= 10)",
"Small (1 < X <= 3)",
"Singlet",
"No TCR information"))
category.order <- c("Hyperexpanded (30 < X)", "Large (10 < X <= 30)",
"Medium (3 < X <= 10)","Small (1 < X <= 3)",
"Singlet", "No TCR information", "No TCR data")
qp <- query$umap %>% as.data.frame %>% cbind(metadata) %>%
ggplot(aes(x=UMAP1, y=UMAP2, color=cloneCategory)) + geom_point(size = 1) +
theme_bw(base_size = 18) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +coord_fixed() +
scale_color_manual(values = colors.no.info)
qp
fig.size(5,8)
metadata %>%
group_by(cell_type_pred_knn, cloneCategory) %>%
dplyr::count() %>%
group_by(cell_type_pred_knn) %>%
mutate(Percent=100*n/sum(n)) %>%
ungroup() %>%
ggplot(aes(x=cell_type_pred_knn, y=Percent, fill=cloneCategory)) +
geom_col() +
# scale_fill_brewer(palette = "Blues", direction = -1, name = "Clone Type") +
scale_fill_manual(values = colors.no.info) +
scale_y_continuous(labels = function(x) paste0(x, "%")) +
ggtitle(NULL) +
labs(fill = "Sample") +
xlab("cell states") + ylab(NULL) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
legend.text = element_text(size=12),
legend.title = element_text(size = 14),
plot.margin = margin(t = 0, r = 0, b = 0, l = 2, unit = "cm"))
metadata %>%
group_by(TYPE, cloneCategory) %>%
dplyr::count() %>%
group_by(TYPE) %>%
mutate(Percent=100*n/sum(n)) %>%
ungroup() %>%
ggplot(aes(x=TYPE, y=Percent, fill=cloneCategory)) +
geom_col() +
# scale_fill_brewer(palette = "Blues", direction = -1, name = "Clone Type") +
scale_fill_manual(values = colors.no.info) +
scale_y_continuous(labels = function(x) paste0(x, "%")) +
ggtitle(NULL) +
labs(fill = "Sample") +
xlab("cell states") + ylab(NULL) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
legend.text = element_text(size=12),
legend.title = element_text(size = 14),
plot.margin = margin(t = 0, r = 0, b = 0, l = 2, unit = "cm"))
# merged@meta.data %>%
# group_by(tissue, cloneCategory) %>%
# dplyr::count() %>%
# group_by(tissue) %>%
# mutate(Percent=100*n/sum(n)) %>%
# ungroup() %>%
# ggplot(aes(x=tissue,y=Percent, fill=cloneCategory)) +
# geom_col() +
# # scale_fill_brewer(palette = "Blues", direction = -1, name = "Clone Type") +
# scale_fill_manual(values = colors.no.info) +
# scale_y_continuous(labels = function(x) paste0(x, "%")) +
# ggtitle(NULL) +
# labs(fill = "Sample") +
# xlab("tissue") + ylab(NULL) +
# theme_classic() +
# theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
# legend.text = element_text(size=12),
# legend.title = element_text(size = 14))
fig.size(8, 12)
#getting clones for each cluster
# unique.states <- unique(merged$cell.states)
unique.tissue <- unique(metadata$tissue)
all.states.clones <- list()
for (s in unique.tissue){
a <- metadata %>% filter(tissue == s) %>%
select(T_clonotype_id) %>% na.omit() %>% distinct(T_clonotype_id)
all.states.clones[s] = a
}
#using the upsetR package
# all.states.clones
UpSetR::upset(UpSetR::fromList(all.states.clones), sets = unique.tissue, keep.order = F,
mainbar.y.label = "Number of Overlapping Clones", text.scale = 2,
sets.x.label = "Unique Clones in tissue", order.by = "freq")
fig.size(8, 18)
#getting clones for each cluster
# unique.states <- unique(merged$cell.states)
unique.tissue <- unique(metadata$cell_type_pred_knn)
all.states.clones <- list()
for (s in unique.tissue){
a <- metadata %>% filter(cell_type_pred_knn == s) %>%
select(T_clonotype_id) %>% na.omit() %>% distinct(T_clonotype_id)
all.states.clones[s] = a
}
#using the upsetR package
# all.states.clones
UpSetR::upset(UpSetR::fromList(all.states.clones), sets = unique.tissue, keep.order = F,
mainbar.y.label = "Number of Overlapping Clones", text.scale = 2,
sets.x.label = "Unique Clones in cell.state", order.by = "freq")
1
fig.size(6,10)
Idents(tregs) <- "cell.states"
DimPlot(tregs, reduction = "umap", cols = cell.state.colors.v2$cell.states, pt.size = 0.8)
fig.size(8, 12)
query$umap %>% as.data.frame %>%
cbind(query$meta_data) %>%
ggplot(aes(x=UMAP1, y=UMAP2, color=cell.states)) + geom_point() +
scale_color_manual(values = cell.state.colors.v2$cell.states) +
theme_bw(base_size = 18) + guides(colour = guide_legend(override.aes = list(size=2))) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + coord_fixed()
fig.size(3,4)
FeaturePlot(tregs, features = "IL2RA", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "FOXP3", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "AREG", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "CXCR6", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "CXCR4", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "CXCL13", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "TCF7", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "NR3C1", reduction = "umap") + coord_fixed()
FeaturePlot(tregs, features = "NR4A2", reduction = "umap") + coord_fixed()
#non expressed - EGF, EPGN, NRG3,
ligands <- c("VEGF", "AREG", "TGFA", "HBEGF", "EREG", "BTC", "NRG4", "NRG1", "NRG2")
ligands2 <- c("NR3C1", "NR4A2", "NR4A1")
fig.size(10,10)
FeaturePlot(tregs, features = ligands, reduction = "umap")
fig.size(3,10)
FeaturePlot(tregs, features = ligands2, reduction = "umap", ncol =3)
Warning message: “The following requested variables were not found: VEGF”
fig.size(5,15)
FeaturePlot(tregs, features = "AREG", reduction = "umap", split.by = "TYPE")
Tregs.markers <- wilcoxauc(tregs, "cell.states")
options(repr.matrix.max.cols=150, repr.matrix.max.rows=200)
Tregs.markers %>%
group_by(group) %>%
filter(padj < 0.05 & auc > 0.6) %>%
arrange(group, desc(logFC)) %>%
#slice_max(n = 40, order_by = logFC)%>%
dplyr::select(feature, group) %>%
mutate(row = row_number()) %>%
tidyr::pivot_wider(names_from = group, values_from = feature) -> Treg.show
Treg.show[1:50,]
| row | AREG+ Treg | CD161+mem. Treg | CD25high Treg | CD25highCXCR6+ Treg | CD25int Treg | GZM Treg | ISG Treg | Naive Treg | Prolif. | TNFAIP3+ Treg |
|---|---|---|---|---|---|---|---|---|---|---|
| <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> |
| 1 | FOS | AQP3 | LGALS1 | TNFRSF18 | RPS29 | GZMK | ISG15 | LEF1 | DEK | CREM |
| 2 | JUNB | GIMAP7 | S100A4 | TNFRSF4 | IFITM1 | CCL5 | MX1 | EEF1B2 | STMN1 | JUNB |
| 3 | DUSP1 | GPR25 | LGALS3 | CTSC | RPS12 | GZMA | IFI6 | TCF7 | TUBA1B | ZFP36 |
| 4 | ZNF331 | CD52 | S100A6 | LINC01943 | LEF1 | IFNG-AS1 | LY6E | RPS3A | BCAP31 | FOSB |
| 5 | JUND | S100A11 | S100A10 | CTLA4 | RPS27 | GIMAP5 | IFIT1 | RPS12 | NME2 | DUSP2 |
| 6 | FOSB | KLRB1 | ANXA2 | BATF | TXNIP | DENR | ISG20 | RPL3 | DUT | RGCC |
| 7 | JUN | CNN2 | EMP3 | CD74 | RPL32 | PRKCQ-AS1 | HERC5 | RPS6 | NDUFA6 | DUSP1 |
| 8 | CREM | FOXP1 | CD74 | TIGIT | RPL30 | GLMN | EPSTI1 | TXNIP | TXN | ZNF331 |
| 9 | ZFP36 | GIMAP4 | HLA-DRB1 | SPOCK2 | RPL13 | TIMM22 | XAF1 | RPL13 | CCDC167 | FOS |
| 10 | DUSP2 | NOSIP | S100A11 | PKM | RPS8 | CHI3L2 | STAT1 | RPS5 | YBX1 | TNFAIP3 |
| 11 | NR4A2 | GIMAP5 | CRIP1 | CXCR6 | RPS18 | ACACA | OAS1 | RPS8 | GAPDH | CXCR4 |
| 12 | CD69 | TRAF3IP3 | ISG20 | IL2RA | RPS21 | SREBF1 | IFITM1 | RPL32 | SUMO3 | DNAJB1 |
| 13 | TNFAIP3 | EMP3 | HLA-DPA1 | CD2 | RPS3A | KRBOX4 | MX2 | RPS13 | PFN1 | LMNA |
| 14 | RGCC | TLE5 | SELPLG | TYMP | RPL34 | NA | SP100 | SELL | C17orf49 | NR4A2 |
| 15 | DNAJB1 | LIMS1 | VIM | GAPDH | RPS13 | NA | IFI16 | CCR7 | PRDX3 | JUND |
| 16 | CXCR4 | SLC25A5 | HLA-DPB1 | ENO1 | RPS28 | NA | OASL | RPL10A | GSTP1 | JUN |
| 17 | SLC2A3 | RPLP0 | LY6E | FOXP3 | RPL8 | NA | SMCHD1 | RPS18 | TWF2 | PHLDA1 |
| 18 | KLF6 | COTL1 | CD52 | TNFRSF1B | RPL10 | NA | IRF7 | RPLP2 | DCTPP1 | DNAJA1 |
| 19 | DNAJA1 | KLF3 | CLIC1 | HLA-DRB1 | RPL3 | NA | TRIM22 | EEF1G | CDT1 | CD69 |
| 20 | BTG2 | ARPC1B | CD99 | LGALS3 | RPL18A | NA | EIF2AK2 | RPL5 | PHF19 | PPP1R15A |
| 21 | YPEL5 | LAPTM5 | SYNE2 | SRGN | RPL11 | NA | RNF213 | RPL9 | OVCA2 | DUSP4 |
| 22 | PPP1R15A | ACTG1 | SH3BGRL3 | CARD16 | RPL39 | NA | SAMD9 | GIMAP7 | ACTB | YPEL5 |
| 23 | RGS1 | MGAT4A | LIME1 | PGAM1 | RPL36 | NA | CMPK2 | RPL29 | PSMA3 | TSC22D3 |
| 24 | SOCS1 | RCSD1 | IL32 | RNF213 | RPS3 | NA | IFI44L | FOXP1 | MDH1 | TENT5C |
| 25 | TENT5C | RPL7 | TLE5 | LGALS1 | RPLP2 | NA | UBE2L6 | RPS23 | TRBV6-1 | SELENOK |
| 26 | ZFP36L2 | IKZF3 | LSP1 | TPI1 | RPL22 | NA | SLFN5 | RPL21 | HLA-DMA | HSPH1 |
| 27 | TSC22D3 | RPSA | CYBA | IL2RB | RPL21 | NA | DRAP1 | RPL4 | SNRPE | KLF6 |
| 28 | FTH1 | MSN | TMSB10 | IL2RG | RPL14 | NA | LGALS9 | RPL7 | RBPJ | BTG2 |
| 29 | UBE2S | RIPOR2 | HLA-A | CST7 | RPL19 | NA | SP110 | RPL22 | MCM5 | SLC2A3 |
| 30 | SELENOK | COX6C | HLA-C | MYL6 | RPS23 | NA | SAMD9L | RPL18A | PSMB2 | RGS1 |
| 31 | NR3C1 | CCR4 | PFN1 | LAG3 | RPS6 | NA | TAP1 | RPS21 | MCM7 | FTH1 |
| 32 | HSPA1B | ABRACL | ARHGDIB | DUSP4 | RPS27A | NA | PARP9 | RPL13A | PDAP1 | NR4A1 |
| 33 | RGS2 | S1PR4 | TMSB4X | PHACTR2 | EEF1B2 | NA | BST2 | EEF1A1 | COA3 | SRGN |
| 34 | HSP90AA1 | CYTH1 | B2M | TIMP1 | RPL41 | NA | IFIT3 | RPL30 | SUCLG1 | SAMSN1 |
| 35 | HSPH1 | RORA | NA | CKLF | RPL12 | NA | OAS3 | RPL36A | PPIA | HSP90AA1 |
| 36 | TAGAP | CCND3 | NA | IGFLR1 | RPL5 | NA | USP18 | RPL34 | MRPL22 | UBE2S |
| 37 | SRGN | RPL4 | NA | PRDM1 | RPL10A | NA | SPATS2L | RPS27A | DDX49 | SOCS1 |
| 38 | TSPYL2 | CDC25B | NA | LAYN | RPL37 | NA | IFI35 | RPL38 | SLBP | MYADM |
| 39 | FAAH2 | CD6 | NA | PTPRJ | RPS15A | NA | IFI44 | RPS2 | CISD1 | BTG1 |
| 40 | CCNH | MKRN1 | NA | S100A4 | RPL29 | NA | OAS2 | GAS5 | RPA3 | CSRNP1 |
| 41 | H3F3B | MVP | NA | DNPH1 | FAU | NA | SYNE2 | RPL39 | CLSPN | NR4A3 |
| 42 | NR4A1 | SYTL1 | NA | ACP5 | RPS14 | NA | ADAR | RPS3 | PTTG1 | NAMPT |
| 43 | PHLDA1 | PLSCR3 | NA | IL12RB2 | RPS5 | NA | TYMP | RPS14 | SLC35E4 | HSPA1B |
| 44 | GEM | ITGB7 | NA | CD4 | RPL38 | NA | RSAD2 | RPL19 | CMC1 | CCNH |
| 45 | TUBB4B | GIMAP1 | NA | CBLB | RPS25 | NA | ZBP1 | RPL36 | POLR3K | PTPN22 |
| 46 | LMNA | HMOX2 | NA | CACYBP | FTL | NA | PSMB8 | LDHB | PCLAF | FOSL2 |
| 47 | MTRNR2L12 | AC004687.1 | NA | SH2D2A | RPL35A | NA | USP15 | RPL37 | APOBEC3G | HSPA8 |
| 48 | DDIT4 | TC2N | NA | TXN | RPL26 | NA | GBP1 | RPS20 | NABP2 | IRF1 |
| 49 | EIF1 | CDC123 | NA | GBP5 | RPL18 | NA | FOXP3 | RPL23A | LIG1 | NFKBIA |
| 50 | SOCS3 | TMSB10 | NA | UCP2 | RPS7 | NA | DDX60L | PABPC1 | MRPS7 | HSP90AB1 |
mg <- read.csv("/data/brennerlab/Shani/projects/Treg/useful_tables/metabolism_genes_manually_curated.csv")
mg %>% head
| Glycolysis | FAO | OXPHOS | |
|---|---|---|---|
| <chr> | <chr> | <chr> | |
| 1 | HK1 | ACADVL | PDHA1 |
| 2 | HK2 | ACADM | DLST |
| 3 | PFKP | HADHA | SDHA |
| 4 | PGAM1 | HADHB | CS |
| 5 | ENO1 | CPT1A | ACO1 |
| 6 | PKM | CPT2 | IDH1/2/3 |
fig.size(4,8)
FeaturePlot(tregs, features = mg$Glycolysis, reduction = "umap") + coord_fixed()
Warning message: “The following requested variables were not found: HK2, PFKP, PGAM1, ENO1, PKM, LDHA, PDHA1, TPI1, GAPDH”
fig.size(8, 12)
FeaturePlot(tregs, features = mg$FAO, reduction = "umap")
Warning message: “The following requested variables were not found: ”
fig.size(6, 10)
FeaturePlot(tregs, features = mg$OXPHOS, reduction = "umap")
Warning message: “The following requested variables were not found: ACO1 , IDH1/2/3, OGDH , SUCLA2/SUCLG1/SUCLG2, FH , MDH1/2, DLD ”
RA <- readRDS("/data/brennerlab/Shani/projects/Treg/analysis/integrated/integrated.Tregs.rds")
length(rownames(RA))
length(rownames(tregs))
inters <- intersect(rownames(RA), rownames(tregs))
length(inters)
tregs@meta.data <- tregs@meta.data %>% mutate(donorID = patient) %>% mutate(tissue = TYPE)
merged <- merge(RA[inters,], tregs[inters,])
merged
An object of class Seurat 30803 features across 24432 samples within 1 assay Active assay: RNA (30803 features, 0 variable features) 2 layers present: counts, data
merged <- NormalizeData(merged, normalization.method = "LogNormalize", scale.factor = 10000, verbose = F) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>%
ScaleData() %>%
RunPCA(verbose = F)
Centering and scaling data matrix
unique(merged$orig.ident)
unique(merged$tissue)
fig.size(4,4)
set.seed(1)
merged <- RunHarmony(merged, group.by.vars=c("orig.ident", "donorID", "tissue"), reduction.use = "pca",
plot_convergence = TRUE, max_iter = 10,
early_stop = T)
Transposing data matrix Initializing state using k-means centroids initialization Harmony 1/10 Harmony 2/10 Harmony 3/10 Harmony converged after 3 iterations
set.seed(1)
merged <- Run_uwot_umap(merged)
merged <- FindClusters(merged, graph.name = 'humap_fgraph', resolution = 0.8, verbose = TRUE)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 24432 Number of edges: 304539 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.7036 Number of communities: 11 Elapsed time: 2 seconds
merged$tissue <- factor(merged$tissue, levels = c("HD.PBMC", "RA.PBMC", "RA.Syn.Fluid", "RA.Syn.Tissue", 'Blood', 'SF', 'Tissue'))
fig.size(8,8)
DimPlot(merged, label = T, repel = T, pt.size = 0.5, label.size = 6) +
scale_color_tableau("Tableau 20") + theme(text = element_text(size = 20))
DimPlot(merged, group.by = "donorID")
DimPlot(merged, group.by = "tissue")
DimPlot(merged, group.by = "orig.ident")
DimPlot(merged, group.by = "cell.states") + scale_color_manual(values = cell.state.colors$cell.states)
fig.size(6, 20)
DimPlot(merged, group.by = "cell.states", split.by = "tissue") + scale_color_manual(values = cell.state.colors$cell.states)
fig.size(8, 12)
tregRef$umap$embedding %>% as.data.frame %>%
cbind(tregRef$meta_data) %>% filter(tissue != "HD.PBMC") %>%
ggplot(aes(x=UMAP1, y=UMAP2, color=cell.states)) + geom_point() +
scale_color_manual(values = cell.state.colors$cell.states) +
theme_bw(base_size = 18) + guides(colour = guide_legend(override.aes = list(size=2))) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + coord_fixed()
fig.size(12,16)
tregRef$umap$embedding %>% as.data.frame %>%
cbind(tregRef$meta_data) %>%
ggplot(aes(x=UMAP1, y=UMAP2, color=cell.states)) + geom_point() +
scale_color_manual(values = cell.state.colors$cell.states) + facet_wrap(~tissue) +
theme_bw(base_size = 18) + guides(colour = guide_legend(override.aes = list(size=2))) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + coord_fixed()